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Abstract. We present here a series of numerical simulations of the development of Kelvin-Helmholtz instability in a rela- 
tivistically hot plasma. The physical parameters in the unperturbed state are chosen to be representative of local conditions 
encountered in Pulsar Wind Nebulae (PWNe), with a main magnetic field perpendicular to a mildly relativistic shear layers. By 
using a numerical code for Relativistic MHD, we investigate the effect of an additional magnetic field component aligned with 
the shear velocity, and we follow the evolution of the instability to the saturation and turbulent regimes. Based on the resulting 
flow structure, we then compute synchrotron maps in order to evaluate the signature of Kelvin-Helmholtz instability on the 
emission and we investigate how the time scale and the amplitude of the synchrotron modulations depend on shear velocity and 
magnetic field. Finally we compare our results to the observed variable features in the Crab Nebula. We show that the Kelvin- 
Helmholtz instability cannot account for the wisps variability, but it might be responsible for the time dependent filamentary 
structure observed in the main torus. 
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1. Introduction 

. Pulsar Wind Nebulae (PWNe) arise from the confinement of 
' the ultra-relativistic wind, powered by the pulsar spin-down en- 
, ergy, by the surrounding supernova ejecta, which are in free 
' expansion in the early phases after the supernova event. The 
prototype of this class of objects is certainly the Crab Nebula. 
PWNe are characterized by continuous non thermal emission 
(synchrotron and inverse Compton), due to the injection of par- 
ticles downstream of the termination shock (the discontinuity 
between the wind and the PWN) and to the presence of a mag- 
netic field, which is basically the residual (toroidal) field car- 
ried by the wind, amplified after the termination shock (Kennel 
& Coroniti .1984 1. 

Recent 2D-axisymmetric Relativistic MHD (RMHD) simu- 
lations of PWNe (Komissarov & Lyubarsky 2004 Del Zanna et 
al. l2004> . aimed at reproducing the main dynamical features of 
the inner structure of PWNe such as equatorial rings and polar 
jets, have shown that, due to the anisotropy in the wind energy 
flux (Bogovalov ll999> . high velocity flow channels develop in- 
side such objects. It has been suggested that the interface be- 
tween such channels can be subject to Kelvin-Helmholtz (KH) 
instability (Begelman 1999 1, and that this might be at the ori- 
gin of the time variability of some of the features present in the 
X-ray maps of the Crab Nebula, namely the wisps and possibly 
the filamentary structures in the main torus (Tanvir et al. il997i 
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Hester et al. '2002V However, at present, temporal and spatial 
resolution in global 2D simulations of PWNe is not sufficient 
to properly investigate the development of such instability and 
to follow its growth up to the non linear regime. 

To verify if KH instability can be at the origin of the ob- 
served time variability it is important to understand if the time 
scale for its growth and the signature on the emitted syn- 
chrotron radiation do agree with observations. Global 2D sim- 
ulations also assume a purely toroidal magnetic field (no mag- 
netic field component is aligned with the flow speed), while it 
is known that even a small component of poloidal field might 
substantially affect the growth and evolution of the pertur- 
bation (Miura 1984 Malagoli et al. [1996l leading to differ- 
ent observational signatures. For example, in a purely toroidal 
case the energy cascades toward large wavelengths, while if 
a poloidal field is present the cascade is toward small wave- 
lengths. Synchrotron maps of PWNe seem in fact to suggest 
the possibility of a small scale turbulent magnetic field in or- 
der to explain features as the jet, whose synchrotron emission 
should be strongly suppressed in the case of a purely toroidal 
magnetic field (Komissarov & Lvub arskv I2004> . 

Let us briefly recall how relativity modifies the develope- 
ment of the instability. In the case of a parallel magnetic field, 
relativistic flows tend to stability both in terms of the growth 
rate and range of the unstable velocity (the lower and upper 
cutoff in the Mach number close together. For mildly relativis- 
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tic flows the lower cutoff shifts to lower values (destabiliza- 
tion at smaller speeds) while for ultrarelativistic flows instabil- 
ity is present only for oblique propagation. In general a rela- 
tivistically hot fluid has slightly higher growth rates, but devi- 
ations are important only very close to the stability threshold. 
In the case of a purely transverse magnetic field there is always 
a range of instability even for parallel propagation. Again the 
upper (and lower cutoff) at first decreases for mildly relativis- 
tic flows and then increases. As in the classic case secondary 
modes are found to develop for oblique propagation. The re- 
gion for instability of secondary modes overlaps with that of 
primary modes. It is also found for very oblique propagation in 
the ultrarelativistic case that the instability region splits and a 
narrow range of stable modes is found. 

Concerning the non linear phases, in general one expects 
that in a relativistically hot plasma both magnetic field and in- 
ternal energy contribute to the inertia. This implies that even if 
P » pv^, the incompressible approximation fails. The conse- 
quence is that the pressure fluctuations due to the formation of 
large eddies, can be of the same order of the background pres- 
sure. Concerning the fluctuations of the velocity components, 
it is found that the amplitude of the velocity perturbation in 
the direction parallel to the interface is smaller than that per- 
pendicular to the interface. This is indeed a consequence of the 
relativistic dynamics. 



Recently Bodo et al. J2004> have shown that it is possible 
to reformulate the KH instability for relativistic unmagnetized 
fluid in an equivalent way to the non-relativistic case, by re- 
defining the Mach number. Given that the equations for trans- 
verse relativistic MHD have the same mathematical structure of 
the one for relativistic HD (Romero et al. .2005j . we can con- 
clude that we do expect a similarity to hold also for transverse 
MHD. 

Here we present the first numerical simulations of the 
growth of KH instability in a RMHD regime, in the simplified 
plane parallel geometry, which can be considered an extension 
of previous works by Miura ( 1984 1 and Malagoli et al. (1996 1, 
with which our simulations share many common points. KH 
instability has been widely studied by means of numerical sim- 
ulations in the non relativistic regime, both magnetized and un- 
magnetized (e.g. Keppens et al. Il999l and references therein). 
The recent development of numerical schemes for relativistic 
fluid dynamics and MHD allow us to finally extend such results 
to the relativistic regime. Recently numerical simulations were 
presented for unmagnetized high Lorentz factor jets (Perucho 
et al. 2004 1. Our simulations mainly focus on the subsonic or 
mildly supersonic regime typical of the flow inside PWNe, and 
take into account the possible presence of a poloidal magnetic 
field. Based on the numerical results for the flow structure we 
then compute synchrotron emission maps associated with the 
development of the instability, to be compared with observa- 
tions of the Crab Nebula. We also evaluate the effect of the 
poloidal component on the polarization properties of the emit- 
ted radiation. 



2. Numerical method and initial conditions 

Simulations were performed using the shock-capturing code 
for RMHD developed by Del Zanna et al. (,2003J . to which 
the reader is referred for a detailed description of the algo- 
rithm. The scheme is particularly simple, since solvers based 
on characteristic waves are avoided in favor of central-type 
component-wise techniques (HLL solver). Density, momen- 
tum, energy, and a passive scalar that traces the mass function 
and allows us to evaluate the mixing, are evolved using a con- 
servative approach, while the magnetic field is evolved via an 
ideal induction equation by following the Upwind Constrained 
Transport method (UCT, Londrillo & Del Zanna 2004 1, which 
guarantees that the solenoidal condition is satisfied up to ma- 
chine accuracy if second order is employed for interpolation. 
We solve the RMHD equations in a 2D planar domain with 
Cartesian coordinates (x, y), where the shear layer is located at 

Given the high Lorentz factor of the pulsar wind (10"^ ^, 
Kennel & Coroniti Il984> . at the termination shock the bulk 
motion energy is converted into disordered (thermal) veloci- 
ties, thus inside PWNe we have a relativistically hot plasma 
plpc^ » 1 (with p thermal pressure and p rest mass density). 
We adopt a value p/pc^ - 20 which is high enough to neglect 
the inertia associated with the rest mass of the fluid (this is 
also a typical value in global simulations of PWNe, where the 
Lorentz factor in the unshocked pulsar wind is ~ 100). We use 
an ideal gas equation of state with adiabatic coefficient F = 4/3, 
appropriate for a relativistically hot plasma. Diff'erent values of 
B, are considered, in order to model different intensities of the 
transverse magnetic field (in the following we will refer to this 
component of the magnetic field as "toroidal", Btoi in analogy 
with 2D global PWNe simulations). In addition to this main 
field, we also investigate cases with a (usually much weaker) 
magnetic field component aligned to the interface, (in the 
following we will refer to the component of the magnetic field 
in the x - y plane as the "poloidal" field, Bpoi, again in anal- 
ogy with the PWN case). Here we will measure the relative 
importance of the magnetic fields by giving the values of the 
non-dimensional toroidal and poloidal magnetization parame- 
ters 
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Since in our configuration B^ is the main component, crtor ^ 1 
means a pressure dominated plasma, cr,oi- » 1 means a mag- 
netically dominated plasma, and the condition of equipartition 
is cTtoi- = 1, while we will always take o-poi <s: 1. 

The shear velocity profile in the initial unperturbed config- 
uration is assumed to be: 



y, = -yotanh(3;/a) 



(2) 



where the value of a has been chosen in order to avoid the 
growth of perturbation on the same scale of the cell size, which 
may occur if a sharp transition is used. In fact, in the limit of 
ff — > oo, the system becomes scale free (Ferrari et al. 1980), and 
given that the growth rate is proportional to the wave number, 
numerical noise can be strongly amplified. The main reason to 
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introduce a shear layer of finite thickness is exactly to prevent 
the growth of such modes. At time f = we also impose a 
transverse monochromatic perturbation V, of the form: 

Vy = Wo sin {k,x) exp [-{yllif] (3) 

where Wo (<k Vo) is the amplitude of the perturbation, is its 
wave number and fi a parameter that describes its fading along 
the y-direction. 

In the jc-direction the size of the computational domain is 
taken D - 2nlkx, and periodic boundary conditions are as- 
sumed, whereas in the y-direction boundaries are placed at 
larger distance, ~ +14D, and 0-th order extrapolation is as- 
sumed in order to minimize possible spurious reflections. We 
use a uniform resolution with 180 zones in the x-direction, 
while in the y-direction a stretched grid is employed with a 
uniform resolution of 360 zones in the range {-D,D]. In all 
our simulations we have chosen a - D/lOO, fi - D/IQ, and 
Wo = 0.01 Vo. The value of a has been chosen to be the small- 
est possible in order not to affect the growth rate. In all simu- 
lations we have considered monochromatic perturbations with 
the longest possible wavelength A - D, for numerical conve- 
nience. From now on we normalize all lengths by this quantity, 
so that for example the x domain is [-1, 1] and kxX - 2nx. The 
chosen time unit is the light travel time Djc, so that all veloci- 
ties will be expressed in units of c. 
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3. Flow dynamics 

We wiU discuss here the dynamical evolution of the fluid struc- 
ture and the magnetic field in our simulations. Given that our 
simulations are basically the relativistic extension of those of 
Malagoli et al. ( 1996 1, our discussion will share many common 
points with that work. In this paper we mainly focus on the sub- 
sonic and marginally supersonic cases. Three different values 
for the shear velocity are considered: Vo - 0.1, typical of the 
flow speed in the outer edge of PWNe; Vq - 0.25, which cor- 
responds to a case with Mach number ■x, \\Y^- 0.4, a typical 
situation encountered in the inner high velocity flow channels. 
Concerning the toroidal magnetic field, we consider cases with 
CTtor - 0.2 (pressure dominated regime), (Xtoi- - 1 (equiparti- 
tion), and (Xtor = 2 (magnetically dominated regime). Different 
values of the poloidal magnetic field B j are then adopted to 
investigate its importance in the stabilization of the interface. 

During the initial phase the evolution of the perturbations 
follows with good accuracy the prediction of the linear theory 
by Ferrari et al. ( 19801, as shown in Fig. Q where the expo- 
nential growth is well reproduced. To allow a more straight- 
forward comparison with the analytic theory, where the growth 
rate is defined in a system a rest with one of the two fluids, 
we consider here the particular case = 0, corresponding to 
their Fig. 1(b) for a relativistically hot plasma'. In this simple 
case the value of the magnetic field that must be used in the 
dispersion relation does not change from the fluid to the per- 
turbation frame, thus allowing an easier comparison. The good 

' In the notation by Ferrari et al. (1980") the corresponding adiabatic 
factor is r = 1, the Mach number is M = yS/c, = VS/?, /? = 2Vo, and 
the plotted growth rates must be corrected by the factors Cj = 1 / Vs 
and To" = 1 - 
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Fig. 1. Evolution of the amplitude of the perturbation AV,, - 
0.5(V'v,maA- - Vy^mi,,) as a function of time (Vo - 0.25; trtor = 0). 
Continuous line: numerical result; dashed line: prediction of 
the analytic theory (Ferrari et al. II980> . From top to bottom: 
cTpoi = 0.01, 0.05, 0.1, 0.2. The growth rates found from 
our numerical simulations are, respectively: 0.21, 0.11, 0.16, 
0.08; to be compared with the expectations of the linear theory: 
0.205, 0.108, 0.165, and 0.075. 
Time is in units of D/c. 



agreement suggests that the value we have chosen for a is small 
enough that the growth rate is not affected. In general we have 
verified that the growth rate tends to be smaller than the predic- 
tion of the vortex sheet approximation for larger values of a, 
and that the discrepancies increase for larger shear layer. Such 
effect was already noticed by Perucho et al. (2004) and can be 
easily understood in terms of stabilization by the shear layer 
itself. In general we find that the growth rates derived from our 
numerical simulations agree with the linear theory (both in the 
case of poloidal and toroidal magnetic field) within a few per- 
cent. 
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In the cases with - 0, corresponding to the assump- 
tion that inside the PWN the magnetic field is purely toroidal, 
it can be shown that transverse RMHD basically behaves as 
hydrodynamics (e.g. Romero et al. 2005). Only pressure-like 
effects survive, while magnetic tension (related to the Alfvenic 
modes) is clearly absent. We find in fact that the development 
of the KH instability and its saturation show only a minor de- 
pendence on the magnetization of the fluid (except the change 
in effective sound speed and compressibility). The system is 
seen to evolve to form a single vortex (not shown here), rather 
stable, on the scale of the computational box. This behavior is 
typical in hydrodynamics, where vorticity evolves to the largest 
possible scales, although here we cannot see a proper inverse 
cascade since we have initialized our system right from the 
start with the largest possible wavelength. A pressure and den- 
sity depression forms rapidly at the center of the vortex. The 
variation of these quantities is a function of the shear veloc- 
ity: for Vo = 0.1, 0.25, and 0.4 the pressure saturates to val- 
ues respectively 10, 50, and 75% lower than in the initial un- 
perturbed state. The transverse velocity Vy grows from Wo to 
value a: +Vo, while the perturbation on the parallel velocity Vx 
is somewhat smaller and ^ O.lVo- The time at which satura- 
tion is reached depends on the value of the shear velocity Vo 
too. In particular, the higher Vo the less it takes to saturate: for 
Vo = 0.1 the saturation time is 9, for Vo = 0.25 it is 5, and 
for Vo = 0.4 goes down to about 3-4. On the other hand, the 
value of the saturation time does not seem to depend on the 
magnetization parameter cr^oT- 

As expected, in the presence of an additional (small) 
poloidal magnetic field component the development of the in- 
stability changes qualitatively and there is an excitation of 
small scale turbulence. As a reference case we will consider 
here the case Vo = 0.25. The only major quantitative dif- 
ferences with the cases with higher and lower shear veloc- 
ity is that the stabilization due to a poloidal magnetic field is 
more efficient the lower the value of Vq. However, even in the 
case Vo = 0.4 a magnetic field such that o-poi - 0.1 is suffi- 
cient to strongly suppress the KH growth. The linear phase of 
the growth in the presence of a weak poloidal magnetic field 
(cpoi = 0.01), is essentially the same as in the hydrodynamical 
case (see Fig.|2j. The magnetic field frozen into the fluid gets 
stretched and twisted as the large scale vortex develops. During 
the linear phase the amplitude of the perturbation in the veloc- 
ity field is approximatively the same as in the hydrodynamic 
(or transverse magnetohydrodynamic) case. Looking at Fig.|3] 
we observe that in the linear phase the poloidal velocity can 
reach values as high as 0.4 (to be compared with Vo = 0.25). 
However, the decay of the MHD instability strongly reduces 
the perturbation and at the center of the shear layer the poloidal 
velocity drops to value close to 0. On the contrary the density 
and pressure perturbations are strongly reduced with respect to 
the hydrodynamical case, to values not exceeding 25%. 

During the linear phase the poloidal magnetic field is am- 
plified to values of order ^j^pVl (equipartition with the kinetic 
shear energy). The amplification does not seem to depend on 
the value of the toroidal magnetic field, probably because in 
our cases the enthalpy is pressure dominated. The other inter- 



esting property is that the poloidal magnetic field perturbation 
reaches a maximum after the time at which the velocity per- 
turbation peaks, with time scales longer of a factor 2, approxi- 
mately. This is probably related to numerical resistivity effects. 
It is know in fact that the higher the resistivity the lower is the 
amplification of the poloidal field. We want also to point out 
that in our 2.5D approximation such magnetic amplification is 
only a transient effect, because in our geometry real dynamo 
processes are not possible. As we can seen from Fig. 12 the 
magnetic field is stretched by vortical flows, the gradient scales 
get smaller, until the resistive scale is reached. This is the mo- 
ment when reconnection of the poloidal magnetic field starts 
and represents the end of the linear phase. This reconnection 
process takes place in an intermittent manner After the first re- 
connection, the magnetic field is again twisted and wound up 
until another reconnection event occurs. Such behavior is also 
reflected in the velocity perturbations. After the first reconnec- 
tion, the velocity perturbation is allowed to grow again. This 
phase is also characterized by an increase in the extent of the 
mixing layer. 

As the magnetic turbulence decays in this direct cascade 
to smaller scales a new statistically steady flow is reached. A 
large mixing layer is formed, whose size depends on the value 
of the shear velocity and is of order of the size of the origi- 
nal perturbation. The mixing layer is characterized by filamen- 
tary structures, which are known features of decaying MHD 
turbulence, approximatively aligned with the direction of the 
unperturbed interface. In comparison with the corresponding 
hydrodynamical case the extent of the shearing layer is about 
twice as bigger Looking at Fig.|2lwe notice that the poloidal 
magnetic field is strongly amplified in those filamentary struc- 
tures (which, as pointed by Malagoli et al. ,1996. turn out to 
be slow magnetosonic waves) to values that, at the end of the 
linear phase, can exceed the amount of toroidal magnetic field. 
However, in the equipartition case cr,or - 1, the poloidal field 
compressed in the filaments tends to saturate, to values from 
10% of the toroidal field for Vo = 0.1, up to 40% for Vo = 0.4 
(apparently with a linear scaling). In the case of stronger field, 
the system tends to evolve to a stable configuration where the 
shear layer is characterized by a series of filamentary structures 
that are slightly inclined with respect to direction of the initial 
interface. 

If we now follow the evolution of the shear layer width, 
we notice that it increases during the linear and subsequent 
reconnection phase, and then saturates at later times. The de- 
gree of mixing inside the layer is complete in the case of a 
weak poloidal field, although there is no appreciable mass loss 
from the layer This implies that if KH instability is at work in 
PWNe, the amplified small scale magnetic fields tend to be con- 
fined in the vicinity of the shear layers. Moreover the extension 
of the mixing layer (evaluated by using a passively advected 
scalar) turns out to be smaller (about one half) than the size of 
the shear layer as defined by the Vv profile, averaged in the x di- 
rection (Fig. l4l5> . With a stronger poloidal field, the growth of 
the shear layer is not accompanied by any mixing. Apparently 
the KH instability is able to amplify the poloidal magnetic field 
on timescales of order of the saturation time. Thus, even if, in 
the simplified geometry of out simulations, permanent dynamo 
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Fig. 2. Evolution of the magnetic structure (Vo - 0.25; crior - 1; (Xpoi = 0.01). The contours represent poloidal magnetic field 
lines while the colors indicate the value of Bpoi/Btor- The wavelength of the initial perturbation is /I = 1. From top to bottom and 
from left to right panels are uniformely spaced in time from 5 to 30. The first panel corresponds to the end of the linear phase. 
The second panel corresponds to the maximum amplification of the poloidal magnetic field before reconnection starts. Notice, at 
later times, the development of the filamentary features discussed in the text. 
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Fig.3. Evolution of the poloidal velocity for the same run (Vo - 0.25; crtor = 1; cTpoi = 0.01) and at the same times as in Fig.|2l 
The maximum in the value of the perturbations in the poloidal velocity is reached at the end of the linear phase, first panel. As 
the system evolves in time, a slow velocity central channel is formed, characterized by a filamentary wave pattern. 
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Fig. 4. Evolution of the integrated velocity profile across the shear layer (Vo - 0.25, crto^ - 1). From left to right: o-poi - 0.01, 
o"poi = 0.05, and cTpoi = 0.1. Note that for f > 15 the profile is relaxed to an almost steady structure, in all cases. For stronger 
poloidal fields the velocity profile evolves to narrower and smoother transitions, typical of the stabilizing effect of the magnetic 
field parallel to the shear layer 



processes are prevented, the decay of the amplified magnetic 
field appears to be very long. 

The amplification of the poloidal magnetic field at shear 
layers may be very important on the observational side as well. 
We recall that, during the linear phase, the field components 
increased by the instability tend to reach equipartition with the 
flow itself, thus for stronger velocity shears a higher poloidal 
field is expected. As pointed out in the introduction, a purely 
toroidal field seems unable to explain emission features like 
the jet. However, material flowing into the jet comes from high 
velocity channels, diverted toward the axis by hoop stresses in a 
sort of back-flow. In the case of a striped wind, that is when we 
have a region with vanishing toroidal field around the equator, 
the shear of such channels is against an equatorial (unmagne- 
tized) flow, which is not inhibited as in the case of full mag- 
netization and that can be actually quite strong (Del Zanna et 
al. ^004). Therefore, the resulting shear layer may contribute 
greatly to the amplification of a small scale poloidal magnetic 
field, which can then in part be advected with the fluid toward 
the axis and thus be responsible for the observed polar emis- 
sion. 

4. Synchrotron emission 

It is evident from the above discussion on the dynamics arising 
from the development of KH instability, that the resulting vor- 
tical motions and the amplification of poloidal magnetic field 
lead to substantial deviations from the laminar unperturbed 
flow. It is natural to expect such modifications to be reflected in 
the emission properties and to provide a possible explanation 
for the observed time variability. Following Bucciantini et al. 
J200 5aV to which the reader is referred for a derivation of the 
equations, the local synchrotron emissivity s(v) at a frequency 
V can be approximated by: 

s(v) oc pD^^Hs'^^^v-^, (4) 

where p is the thermal pressure, 2^ H- 1 is the power law index 
of the particle distribution (in the following we have adopted 



the value ^ = 0.6), B'^ is the magnetic field component perpen- 
dicular to the line of sight, as measured in the fluid frame, and 
D - {y{l - P ■ «)] Here )8 is the flow velocity and n is the 
versor pointing toward the observer Given the magnetic field 
structure one can also compute the polarization fraction and the 
polarization angle. The value of the Stokes parameter U and V 
is derived following Lyutikov et al. ( 2003 ) (see also Bucciantini 
et al. 2005b Blandford & Konigl 1979). To correctly evaluate 
the emission properties it is fundamental that the flow structure 
is computed using relativistic fluid dynamics, and the Doppler 
term D may become one of the most important ingredients. 

Looking at Eq.|4]it is evident that there are two main source 
for modulation of the synchrotron emission. The first is the 
compression or rarefaction connected with the development of 
vortexes, which reflects itself in the enhancement or decrement 
of the p and terms. In the cases with a purely toroidal mag- 
netic field the emissivity in the center of the largest scale vor- 
tex can be strongly suppressed with respect to the surrounding 
laminar flow. In the case Vo - 0.4, for example, the central 
emissivity is about 10 times lower than in the unperturbed re- 
gions. However, in the presence of a poloidal field, the vorticity 
decays to smaller scales (the large scale eddy that forms during 
the linear phase is disrupted, and the resulting turbulent shear 
layer contains many eddies with sizes smaller thn the original 
perturbation) and such compressible modulation is strongly re- 
duced. The second source of modulation is due to the boosting 
term D, which amplifies the emission from a turbulent fluid at 
points where the local velocity is toward the observer Even in 
the presence of a weak poloidal field, the transverse velocity 
Vy can reach values comparable with Vo- The boosting term 
can account for fluctuations in the emissivity of about a factor 
10 for Vo - 0.4 and a factor 4 for Vo - 0.25. Finally, as we can 
see from Fig.|2 the poloidal magnetic field can be amplified to 
a non negligible fraction of the toroidal field. A poloidal field 
on small scales, while still contributing to the emissivity, might 
lead to a depolarization of the observed radiation (if the field is 
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Fig. 6. Time evolution of the synchrotron emission, normalized to the mean value, for an inclination angle of 30°. A purely 
toroidal magnetic field at equipartition crjor = 1 is assumed in the unperturbed state. From left to right: Vq = 0.1, Vo = 0.25 and 
Vo = 0.4. 
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Fig. 5. Evolution of the integrated velocity profile across the 
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panel Vo = 0.4. Notice that for Vo =0.1 the extent of the shear 
layer is smaller than in the case Vo - 0.25 (Fig. |^, while no 
major differences are shown in the Vo - 0.4 case. 



purely toroidal the polarization is by definition at its theoretical 
maximum). 



However, in spite of the local fluctuations in the emissivity, 
associated with the growth of the KH instability, which could in 
principle account for large modulation in the synchrotron emis- 
sion, it is not obvious how, once the contribution from dififerent 
fluid elements is integrated along the line of sight, the result- 
ing modulation will actually look like. Moreover, in the case of 
the Crab Nebula, the inclination of the equatorial plane (which 
should be parallel to the shear interface of the flow channels) 
and the observer is about 30° (Weisskopf et al 2000 1. Once the 
contribution of the various fluid elements along the line of sight 
is integrated, the combined eff'ect might lead to a suppression 
of the synchrotron modulation with respect to what could be 
inferred by just considering local conditions. To compute the 
total emission, the integration along the line of sight was per- 
formed on a layer extended in the range - 1 < y < 1 . In general, 
larger integration layers yield smaller modulations. Such size 
was chosen to guarantee that in all cases the shear layer was 
contained within the integration region. For example, by con- 
sidering just the contribution of the compressible term in Eq.|3 
in the case of a purely toroidal field we find that, for Vo - 0.4, 
the fluctuations in the integrated emission are about 50% for an 
inclination angle of 90°, but they drop to less than 5% for an 
angle of 30°. 

Let us first consider the case of purely toroidal magnetic 
field. Looking at Fig.|6]we observe that the maximum in the 
synchrotron modulation coincides with the end of the linear 
phase (f - 5). Notice also that at this time the synchrotron 
perturbation has a wavelength which is about one half that of 
the original perturbation. This is due to the fact that the com- 
pressible and boosting terms do not sum up in phase. The am- 
plitude of the modulation in the synchrotron emission is how- 
ever extremely small: it increases for higher Vo but in all our 
cases is below 10%. At later times, in the saturation regime, the 
synchrotron modulation is almost completely due to the com- 
pressible terms. Once a stable vortex is established, the effects 
of boosting from different fluid elements tend to compensate 
each other We verify that such a behavior does not depend on 
the assumed inclination angle. These results also show that the 
timescale for the onset of synchrotron modulation are of order 
of the saturation time. 

If a poloidal magnetic field is present the decay of the turbu- 
lence to small scales qualitatively changes the emission signa- 
ture. Until the end of the linear phase and before reconnection 
takes place the modulation in the synchrotron emission is the 



N. Bucciantini and L. Del Zanna: Local Kelvin-Helmholtz instability and synchrotron modulation in Pulsar Wind Nebulae 



9 



same as in the cases with purely toroidal magnetic field. After 
reconnection, the decay to small scales leads to a more turbu- 
lent and time dependent behavior. Looking at Fig. we see that 
the amplitude of the modulation is larger, even if still less than 
10-15%. There is also an intrinsic time variability on timescale 
of a few units (D/c). Modulations do not appear stationary in 
the original perturbation frame: looking at emission peaks, they 
appears to move with a speed ~ 0.1c, smaller than the shear ve- 
locity Vo- This is an important diff'erence. In the purely toroidal 
case (the modulation is stationary in the perturbation frame) 
the motion of the observed features can in principle be related 
to the underlying flow structure, while if a poloidal field is 
present such correlation is not possible. For higher values of 
the poloidal magnetic field the amplitude of the modulation at 
tiie end of the liner phase is higher but, as expected, the sta- 
bilization of the layer leads to a suppression at later times. In 
general, the modulation increases in the linear phase and is later 
dumped. Contrary to the purely toroidal case, where the modu- 
lation in the emission was mainly due to the compressible term, 
now pressure fluctuations are much smaller and the modula- 
tion in the synchrotron maps are mainly due to the boosting 
term. We also computed the polarization properties of the ra- 
diation, to verify if the amplified poloidal field could possibly 
lead to any depolarization. Despite the high value reached by 
Bpoi/Btor, no depolarization was observed. The reason can be 
understood by looking at Fig.|2l the poloidal field is amplified 
to high values only in thin filamentary structures. The contri- 
bution of these regions to the integrated emission is negligible, 
and the resulting polarization is the same as in the unperturbed 
state. 

Let now apply the above results to the case of the Crab 
Nebula. Recent observations have shown intrinsic variability 
of the wisps region and X-ray images have shown the pres- 
ence of time dependent filamentary structure in the main torus 
(the wavelength of the filaments is about 0.2 Ly) (Hester et al. 
l2002t . Recent papers published on the subject were more fo- 
cused on the time variability of the wisps region, and no quan- 
titative information about the observed synchrotron modulation 
is available, either for the wisps or for the torus. It is however 
evident from the maps that wisps are characterized by large am- 
plitude fluctuations (~ 100% of the wisps intensity) while the 
filamentary structure in the main torus is much fainter. If one 
considers the amplitude of the modulation, it is evident from 
our results that fluctuations lager than 10-20% can be hardly 
due to KH instability. We thus conclude that KH is unlikely to 
be at the origin of the variability in the wisps, while it is still 
possible that it could account for the filamentary features ob- 
served in the main torus. 

Moreover, by looking at the time scale for the growth of the 
modulation, we find that this is too long to explain the wisps re- 
gion variability (which in recent PWN models turns out to be 
at the termination shock). Perturbations having wavelengths of 
order of the main torus filaments are expected to grow with 
timescales of about 1-2 yr, in the linear phase. Assuming a 
flow velocity of order c/3 in the equatorial plane (derived from 
Doppler boosting of the torus), perturbations will have moved 
to a distance of about 3-6 ly, which is in rough accord with the 
location of the main torus (this also implies that the amplitude 
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Fig. 7. Time evolution of the synchrotron emission, normal- 
ized to the mean value, for an inclination angle of 30°, a shear 
velocity Vo = 0.25, a toroidal magnetic field at equipartition 
(Ttoi - 1. Upper panel: poloidal field o-poi - 0.01. Lower panel: 
o-poi - 0.05. 



of the original perturbation cannot be too small). Actual data 
do not allow to understand if the observed time variability of 
modulation can be ascribed to a poloidal field, however the in- 
trinsic time variability in our simulations is of order of years, 
much longer than the observed week-month timescales. It is 
more likely that the time variability to be due to the propaga- 
tion of the perturbation, instead of its decay to smaller scales. 

The last point concerns the polarization. Optical data 
(Schmidt et al. 11979 ) show that polarization in the wisps and 
main torus is about 50-70% of the theorical maximum. Such 
depolarization could actually be due just to foreground and 
background emission. However, if we assume depolarization 
to be intrinsic, our result then shows that small scale magnetic 
fields must be injected at the shock, because amplification of 
poloidal field by the KH instability does not produce any ap- 
preciable depolarization. We conclude that KH instability is 
unlikely to be the cause of the wisps variability (it is more 
likely due to global modulation in the PWN flow structure, 
Komissarov & Lyubarsky 2004 Bogovalov et al. 2005 1, while 
it remains a good candidate for the origin of the time dependent 
filamentary structure in the main torus. 
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5. Conclusions 

In this paper we have extended previous results about KH in- 
stability in the presence of a component of the magnetic field 
pai-allel to the interface (Miura fT^MI Malagoli et al. ll996> . to 
a relativistic regime typical of the conditions inside PWNe. For 
the first time, based on our numerical models we have derived 
the eff'ect of KH instability on the emitted synchrotron radia- 
tion and compared the results with the observed time variable 
features in the wisps region and in the main torus of the Crab 
Nebula. 

Regarding the evolution of the instability, we recover a 
similar behavior as in the non relativistic case. For a purely 
toroidal field the vorticity decays to larger scales (the system 
is seen to evolve to form a single vortex rater stabel on the 
scale of the computational box), transverse velocity increases 
to a value close to the shear velocity and pressure and mag- 
netic field drop substantially in the vortex. The saturation time 
is longer for smaller shear velocity, and in the case of our initial 
perturbation is of order of a few D/c.lfa small poloidal field is 
also present, the evolution of the instability changes. The linear 
phase evolves as in the case of a purely toroidal magnetic field, 
until the twisted poloidal component experiences reconnection 
events, which take place in an intermittent manner MHD tur- 
bulence develops and decays to the dissipative small scales set 
by viscosity and magnetic resistivity (in our case the numeric 
viscosity and resistivity of the code), these are the same scales 
on which reconnection takes place. Fluctuations in the pres- 
sure become negligible at this stage, while those in the trans- 
verse velocity are still high. The poloidal magnetic field is com- 
pressed in filamentary features up to values about 10-40% of 
the equipartition field. There is also an increase in the extent of 
the turbulent shear layer 

As we have shown, fluctuations in the local emissivity due 
to KH instability can be very large, however once the emis- 
sion is integrated along the line of sight the eff'ects of different 
fluid patches tend to compensate, and only minor fluctuations 
remain. This stresses the importance of correctly computing 
the contribution from diff'erent fluid elements for the integrated 
emission, when comparing models to observations. Even if lo- 
cal variation can explain the observed features, it is not obvious 
how this will affect the total observed emission. We find, in our 
settings, that in general fluctuations do not exceed 10-20%. In 
the purely toroidal case these are mainly due to compression 
effects on the thermal and magnetic pressure terms, while, in 
the presence of an initial poloidal field, modulations seem to 
be more related to the boosting term, which may turn on due to 
the turbulent velocity field produced after the saturation phase. 
The time scale for the growth of the synchrotron modulation 
is of the same order of the time scale for the instability to sat- 
urate. In the turbulent phase the time scale is still of order of 
few D/c. By comparing our results to observations of the Crab 
Nebula we concluded that local KH instability is unlikely to be 
at the origin of the wisps variability, while it remains a possible 
explanation for the time dependent filamentary structure seen 
in the main torus. 
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